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Abstract.  In  1993,  Li  and  Mayo  gave  a  finite-difference  method  with  second  order  ac¬ 
curacy  for  solving  the  heat  equations  involving  interfaces  with  constant  coefficients  and 
discontinuous  sources  [Proc.  Symp.  Appl.  Math.  Vol.  48,  W.Gautschi  ed.,  AMS,  1993, 
p  311-315].  In  this  paper,  we  improve  the  above  result  by  presenting  a  finite-difference 
method  which  allows  each  coefficient  to  be  taken  different  values  in  different  subregions 
divided  by  the  interface,  that  is  useful  in  applications.  Our  method  also  has  second  order 
accuracy. 


1.  Introduction 

Consider  the  heat  equation 


ut  {P‘U'x)x  +  {Puy)yi  (1.1) 

where  t  £  [a,  oo)  and  (m,  y)  £  fl,  a  rectangular  region  in  R?  with  an  irregular  interface 
T  which  divides  O  to  two  subregions  Qi  and  02.  The  solution  u(x,y,t)  and  its  normal 
derivative  un(x,y,t)  crossing  the  curve  T  are  known  to  be  discontinuous: 

N  =  u+  ~  u~  =  u(x(s),y(s),t ),  (1.2) 

K3  =  ut~uZ  =  9(x(s),  y(s),t),  (1.3) 

where  s  is  a  parameter  of  T,  the  superscripts  +  and  -  denote  the  limiting  values  of  a 
function  from  one  side  in  Q+  and  another  side  in  Q-  respectively. 

In  1993,  Li  and  Mayo[3]  gave  a  finite-difference  method  with  second  accuracy  for 
solving  (1.1),  assuming  that  /  is  discontinuous  but  j3  is  constant.  In  this  paper,  we  present 
a  finite-difference  method  for  solving  (1.1),  which  allows  /3  to  be  taken  different  values  in 
different  subregions  divided  by  the  interface,  and  which  is  useful  in  applications.  More 
precisely,  we  assume: 
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(1,4) 


(x,y)  G  f2+, 
y)  s  fi  , 


where  P+,{3  are  constants  which  can  be  distinct. 

We  organize  the  paper  as  follows:  In  Section  2,  we  establish  local  coordinate  systems 
around  the  interface.  In  Section  3,  we  give  the  correction  terms  for  the  finite-difference 
method.  In  Section  4,  we  show  that  our  method  is  second  accurate.  Finally,  in  Section  5, 
we  give  some  numerical  examples,  in  which  the  actual  solutions  are  known,  to  confirm  the 
theoretical  result. 
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2.  Local  Coordinate  Systems 

We  first  give  local  coordinate  systems  along  the  interface  V  as  [2,3].  When  a  point  on 
the  interface  is  fixed  for  the  origan  ,  we  use  the  normal  direction  as  the  ^-direction,  which 
has  an  angle  $  with  the  x-axis.  Rotating  the  ^-direction  by  ninety  degrees  anti-clockwise, 
we  obtain  the  77-direction.  Now  we  express  the  curve  T  as  a  function  of  the  independent 
variable  77  locally: 


£  =  x(v)- 

(2.1) 

We  express  (1.2), (1.3)  locally  by  using  the  77  coordinate: 

[u]  =  u+  —  u~  =  u(r},  t), 

(2.2) 

[un]  =  u+-u~  =g(r},t), 

(2.3) 

where  u,  g  are  known  in  advance.  From 

s-M- 

(2.4) 

J  i  +  x? 

Using  differentiation  w.r.t.  77,  noting  y(0)  =  0,  we  have: 


K]  =  «n»  K]  =  9,  (2.5, 2.6) 

[Uvv)  =  ~9Xt)T]  +  OJjjij,  [w^Tj]  =  &T]XvV  +  9v 


2 


(2.7, 2.8) 


*  Changing  (1.1)  locally,  we  have:' 


“l  =  +  Pum  +  /(£,  i),  t ), 


(2.9) 


which  implies 


where 


[ucd  —  Kdi  + 


'ut 

.0 


]  1  9Xr)r]  &t]ti 


(2.10) 


(2.H) 


In  (2.5)-(2.11),  all  functions  in  the  right-hand  sides  are  known  except  [^]  in  (2.10)  which 
will  be  explored  in  the  next  section.  p 


3.  Correction  Terms 

We  first  discretize  in  both  of  the  ^-direction  and  the  y-direetion  by  mesh  size  h: 

1  , 

Uxx  ~  OxUitj  —  ^2  (Ut-lj  -  2 Ui'j  +  ui+ij),  (3.1) 

Uyy  «  SyUiJ  =  —(UiJ-l  -  2 Ujj  +  UjJ  +  l).  (3.2) 

We  group  all  the  grid  points  in  Q  into  two  sets.  The  set  Sreg  consists  the  regular  points, 
each  point  in  one  subregion  has  no  neighbor  point  in  the  another  subregion.  The  set  Sirr 
consists  the  irregular  points,  each  point  in  one  subregion  has  at  least  one  neighbor  point 
in  the  other  subregion.  For  a  regular  grid  point,  the  local  truncation  error  of  (3.1), (3.2) 
from  fiuxx  4"  0yy  +  /  is  OQi2').  For  an  irregular  grid  point,  we  need  add  some  correction 
terms  in  (3.1), (3.2)  such  that  the  local  truncation  error  is  0(h),  therefore  the  global  error 
of  the  solution  for  solving  the  heat  equation  is  0(h2)  after  the  discretization  of  time  t  in 
certain  way. 

At  first,  we  relate  the  jumps  w.r.t,  x  and  y  to  the  jumps  w.r.t.  |  and  rj: 

[«*]  =  [u?]cos0  -  [. un)sin6 ,  [%]  =  [u(]sin9  +  [uv]cos9,  (3.3, 3.4) 


cos26,  (3.5) 

where 

\uxx\i  —  [u^]jcos20  —  2[u^r}]cos9sin9  4-  [ur]7j]sin29, 

sin29 ,  (3.6) 

where 

K/di  =  [ti£dis*n20  +  2  [u£n]cos9sin9  +  [u7]T)]cos29. 


[Uyy]  ~  Uyy  1  + 


Ut 

[0 


Uxx,  —  [ UXx ,  1  + 


Ut 

0 
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By  (2,5)-(2.11),  all  the  terms  in  the  right-hand  sides  of  (3.3)-(3.6)are  known  except  [at]. 
Now  we  consider  an  irregular  grid  point  in  ^-direction,  there  are  four  cases:  P 

(a)  ( Xi ,  yj)  €  O”,  (xi+1,  yj)  €  0+;  ( b )  (x *,  Pj)  €  0+,  (xt-x,  pj)  6  0~; 

(c)  (xi,pj)  €  0+,  (xi+i,pj)  6  O  ;  (d)  (xi,pj)  6  0~,  (xi-i,yj)  e  0+. 

For  case  (a),  let  the  intersection  of  the  line  segment  connecting  (xi,pj),(xi+1,yj)  and  T  is 
(x*,pj).  Using  Taylor’s  series  around  x*,  we  have: 

^2  Vj)  -  2u(xh  pj)  +  u(xi+ 1,  g/j-)} 

=  p{[«]  +  M(xi+ 1  -  ar*)  +  ^l(ari+i  -  a:*)2}  +  i£x  +  0(h), 

which  implies 


$xUitj  "H  C x 


[^‘](a'i+l  X*Y 


+  0(h), 


where 


CxUij  —  {[u]  +  [«a;](a:i+i  —  x*)  +  -  (xt+i  —  £*)2}- 


Similarly,  for  case(b),  we  have 


where 


uXx  —  SxUij  +  Cxiiij  H  —  ■— — - b  0(h), 


CxUi,j  =  -  x*)  +  -  x*)2}. 


For  case(e),  we  have 


where 


Uxx  —  $xui,j  ~b  CxUi;j  -)  — 2 - h  0(h), 


CzUi'j  =  ±{[u]  +  [ti*](*i+1  -  *•)  +  -  a:*)2}. 


For  case(d),  we  have 


[f  Kari-x  -  a:*)2 

W®®  —  $x ^i,j  "b  CxUij  ■— j - 1-  0(h), 


(3.10) 


where 


CxUij  =  -£{[«.]  +  [uJ(Xj_i  -  X*)  +  hgli(x(_,  -  X*)2}. 
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Analogously,  in  ^-direction,*  we  also  have  four  cases.  We  add  correction  terms  such 
that  the  local  truncation  is  0(h). 

Now  using  these  correction  terms  in  both  x-  and  y-  directions,  we  obtain  a  system  of 
ordinary  differential  equations: 

Ki)*  =  P{8«Uij  +  SyVij}  +  Si, j  +  P  ]T  {cxuij  +  [  ^  j  — — °~g2*)2c0g2g} 


{°yui,3  +  ~o 

{x*  iVJq  )SSirr 


(xi0  }V*)€Sirr 

M  ry;,  (yjri  -  y*)2sin20 


}  +  0(h), 


(3.11) 


where  i0  =  i  -  1  or  i  +  1,  ja  =  j  -  1  or  j  + 1.  txo,  r„0  =  1  or  -1  according  to  (3.7)-(3.10). 
We  have 

_  f  1 1  ,  [ut]  .  f  1 1  [ut\ 

J  J  =  [^j  +  JT  “>  b  +  jr,  (3.12,3.13) 


which  imply 


—  (ui,j)t  -p  +  0(h). 


(3.14) 


Using  (3.14),  we  have  the  following  system  of  ordinary  differential  equations: 

~  F(ui,j—h  ui,ji  ui+l,j ,  ui,j+l)  = 

+  SyUjj  +  E(g<0  ,r)gs<rr  CxUi,j  +  E(x»,yJ0)eSirr  CyVjj)  +  fi,j 

1  _  V1  .  Qrn^n^-o-**)20032^  «r 1 1 twjo (Wo -y*)2sin2e‘  (3-15) 

^(xi0,ym)€SirrPlj3l  2 ft*  ^(xm,yj0)eSirrP\-pi  - 2 P - 

At  a  regular  grid  point,  the  local  truncation  error  of  the  right-hand  side  of  (3.15)  from 
fiv>xx  fiyy  4"  /  is  0(h2).  In  the  next  section,  we  will  show  that  at  an  irregular  grid  point, 
the  local  truncation  error  is  0(h). 

Finally,  we  discretize  time  t  by  choosing  At  =  h2.  We  use  Crank-Nieholson  method: 

Ui,j,k+i  =  +  At(0.5F(«*  j  _  1  (fc ,  «i_  Uij)k,  ui+ij>k,  ui>j+1>k) 

+0.5F(UjJ-_iifc+i,  u»— 1  ,j,k+li  ui,j,k+li  ui+l,j,k+lj  ui,j+l,k+l))i  (3.16) 

which  implies  the  local  truncation  error  for  discretizing  t  is  0((At)2).(see  [1,4].)  To  solve 
ui,j,k+ 1  from  (3.16),  we  use  S.O.R.  iteration  with  certain  parameter  u>: 


Ui,j,k+ 1  —  ui,3,k, 


(3.17) 


Ui,j,k+ 1  (1  W)'Ui,j>+l  U(Uh3,k  +  0.5F(tijj_  i>k,  ui,j,ki  ^ij+l,fc) 

+  0  5 F(t^n.+1*  u{n+1)  «(**)  »,(n)  ,,(n)  \\  „  _  i  9 

-f-  U.O/*  Ui,j,k+V  Ui+l,j,k+V  ti*,i+1,fc+l))l  71  -  ll2! 

(3.18) 
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’4.  Accuracy  Analysis 

We  first  show  that  the  denominator  of  the  right-hand  side  of  (3.15)  is  bounded  below 
and  above  by  positive  constants: 

Lemma  4.1.  Let  Dij  = 

1  -  Y"  ft  |T|  Txi°  ^Xi°  ~  x*}2cos2e  _  sr  «  F1!  %o^io  “  V*)2sin26 


Then 


min(0+,0~)  fl 

mJ(/3+,, S-)  -  D'J  -  l  +  max(fl+,fl~)  i 


where  i0  -  i  - 1  or  i  + 1,  j0  -  j  - 1  orj  +  1.  tx.q,  Ty.Q  -  1  or  -1  according  to  (3.7)-(3.10). 

Proof:  At  first  we  prove  the  lower  bound.  Look  at  the  first  summation.  In  x-direction, 

we  have  four  cases. 

(a)  (xuVj)  <=  Q“,  (xi+1,yj)  €  0+.  Then  xio  =  ari+1,  rXiQ  =  -1,0  =  0~, 

and  0  rx.Q  =  1  —  which  is  positive  iff  0~  <  0+. 

(b)  (xhVj)  e  Q+,  (®<_i,  jy)  €  fi~.  Then  xio  =  Xi_ls  rXi(}  =  1,  0  =  /3+, 
and  /3  rXio  =  1  -  |^  which  is  positive  iff  0+  <  0~. 

(c)  (®t,2/i)  €  Q+,  (xi+1,yj)  €  Then  xio  =  xi+1,  rXiQ  =1,0  =  0+, 

and  0  rXiQ  =  1  -  which  is  positive  iff  0+  <  0~. 

(d)  (xuyj)  €  Q~,  (xi+1,yj)  e  0+.  Then  xio  =  x,  rIio  =-1,0  =  0~, 
and  0  ^  rXio  =  1  -  jp-  which  is  positive  iff  0~  <  /3+. 

For  all  cases,  a  term  in  the  first  summation  is  positive  iff  0  =  min{0+,  0~}.  Similarly 
we  can  show  that  a  term  in  the  second  summation  is  positive  iff  0  =  min{0+,0~}  too. 
Only  the  positive  terms  will  reduce  the  lower  bound  of  the  left-hand  side  of  (4.1).  So  the 
left-hand  side  is  bounded  below  by 


Dij  >  1  -  0.5(1  -  -  0.5(1  -  =  rnin(l3+,l3-) 

max(0+,0-)'  K  max(0+,0-))  max(0+,0-)' 

Now  we  turn  to  the  upper  bound  which  can  be  proved  directly  from  (4.1): 

Di,j  <1  +  0.50  —  +0.50  ^  <  1  +  max(0+ , 0~)  i  . 

The  lower  bound  in  (4.2)  is  useful  for  the  stability  of  the  numerical  scheme  and  the 
upper  bound  will  be  used  for  the  following  theorem: 
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'  Theorem  4.2. 


F(ui,j-u  ui-i,j,  Ui,j ,  Ui+itj,Uitj+i)  —  (f3ux x  +  fiktyy  +  f)(xi,  yj)  =  0(h),  (4.3) 

where  (xi}  yj)  is  an  irregular  grid  point  in  Cl. 

Proof:  From  (3.15),  (3.14),  (4.2)  and  (3.11),  we  have: 

Uj  j_|_i)  =  (lijj)t '  = 

f3(6xUitj  +  SyUiJ  +  CXUiJ  +  ^2  CyUiJ )  +  /jj 

(zigiy*)€Sirr  (®*  jS/jq  )€§irr 


-<«*>.  £  e  :r 

\(xi0>ym)£Sirr  (.X*  ,yjo)eSirr  ^  J 

=W»“M  +  SyUij}  +  /„  +  0  E  {<*<«  +  fe]  ~f : >2c°s29} 

+*  E  {W+[|]^.-;;>We}+ow 


(s*,%)6Sirr  L  J 

+  f)(xi,  yj)  +  0(h) 


which  proves  (4.3). 

At  a  regular  grid  point,  the  local  truncation  error  for  discretization  is  0(h2).  At  an 
irregular  grid  point,  the  local  truncation  error  is  0(h)  by  Theorem  4.2.  The  discretization 
of  time  is  0((At)2)  =  0(h 4).  All  these  imply  that  the  numerical  solution  has  global  error 
0(h2). 


5,  Numerical  Examples 

We  choose  some  examples,  in  which  the  actual  solutions  are  known,  therefore  numer-- 
ical  error  computations  can  be  obtained  to  confirm  the  theoretical  result  of  our  method. 
We  choose 

»  0  =  [-1,1]  x  [-1,1],  t  €  [0,  oo],  (5.1, 5.2) 

0+  =  {(a;,  y)  G  Cl  \x2  +  y2  >  1/4},  CT  =  {(a;,  y)  €  Q  |rzr2  +  y2  <  1/4},  (5.3, 5.4) 

r  =  {(as,  y)  6  Cl  \x2  +  y2  =  1/4}.  (5.5) 
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(5.6) 


'The  actual  solution  is  known  in  the  case  /  =  0: 

t  *  1 

u{x,y,t)^-e 

v 

We  give  the  initial  condition  when  t  =  1  and  boundary  condition  when  x  or  y  =  1  or  —  1 
by  using  formula  (5.6).  We  choose  w=1.75  in  (3.18).  For  different  pairs  of  j3+  and  /3~,  in 
t  from  1,0  to  1.5,  we  obtain  the  following  tables,  in  each  the  error  is  computed  by  using 
Euclidean  norm:  ° 


h 

r 

error 

ratio 

0.100 

1000 

i 

2.227782 D  -  04 

- - 

0.050 

1000 

i 

5.391984D  -  05 

4.13 

0.025 

1000 

i 

1.307318D  -  05 

Table  5.1 

4.12 

h 

fi¬ 

error 

ratio 

0.100 

1 

lm 

2.392997 D  -  05 

- -  , — 

0.050 

1 

1000 

6.703319D  -  06 

3.57 

0.025 

1 

1000 

1.766744T>  -  06 

Table  5.2 

3.79 

h 

fir 

error 

ratio 

0.100 

5 

1 

2.629529D  -  04 

- — 

0.050 

5 

1 

6.351060D  -  05 

4.14 

0.025 

5 

1 

1.550294D  -  05 

Table  5.3 

4.10 

h 

fir 

error 

ratio 

0.100 

1 

5 

5.461059D  -  05 

- - 

0.050 

1 

5 

1.309861 D  -  05 

4.17 

0.025 

1 

5 

3.263757 D  -  06 

4.01 

Table  5.4 


The  tables  show  that  when  h  is  reduced  to  a  half  of  it,  the  error  is  reduced  approxi¬ 
mately  to  a  quarter  of  it.  That  confirms  the  numerical  solution  has  second  order  accuracy. 
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